% run this file to recreate Table 2 ("Model Estimation - Targeted and Simulated Moments")
% and Table 3 ("Model Estimation - Estimated Parameter Values")

% Note that this requires repeated model simulations under different values of parameters for the routine trying to match targeted and simulated moments;
% as such the progra will take a while to run (12-24 hours on a laptop)

% The five targeted moments are different whether the *early 2000s*
% case, the *late 2000s* case or the *mid2010s* case are used
% These targets are in the 3,5 matrix called moments_mat
% and the program selets them depending on the case chosen by the user

% Note also that this program was run several times, with several
% alternative initial guesses, to assess sensitivity. The results contained
% in Table 2 and 3 reflect our best assessment of the global optimum
% The program below is a subset of all this work : it shows that using the
% published values as guesses leads to no change in the final values and no
% changes in the distance between targeted and simulated moments



% Define the matrix containing the targeted moments
moments_mat = zeros(3,5);

% moment 1 : fraction of simulated volatility in spot prices attributable to permanent shocks
moments_mat(1,1)= 0.1;
moments_mat(2,1) = 0.5;
moments_mat(3,1) = 0.3;

% moment 2 : s.d. of (HP-filtered) oil inventories, relative to s.d. of (HP-filtered)  GDP
moments_mat(1,2) = 6.0;
moments_mat(2,2) = 5.0;
moments_mat(3,2) = 4.0;


% moment 3 : autocorrelation of GDP
moments_mat(1,3) = 0.815;
moments_mat(2,3) = 0.85;
moments_mat(3,3) = 0.874;

% moment 4 : corr between (HP-filtered) spot price of oil and (HP-filtered) oil inventories
moments_mat(1,4) = -0.4;
moments_mat(2,4) = -0.45;
moments_mat(3,4) = -0.5;

% moment 4 : corr between (HP-filtered) spot price of oil and (HP-filtered)  GDP
moments_mat(1,5) = -0.07;
moments_mat(2,5) = 0.1;
moments_mat(3,5) = 0.36;


poids = eye(5); %how to weigh success in fitting each of the moments

% ready to start estimating

% which case to estimate?
answcase = 3; % select the case to be estimated
              % 1 for *early 2000s*
              % 2 for *late 2000s*
              % 3 for *mid2010s*

% this will be our target
moments = moments_mat(answcase,:)';


% Defining starting values for parameter vector (guesses) for each case

if answcase == 1

thetaguess = [ 0.0655
               0.8158
               0.3694
               0.2707
               0.8466];
elseif answcase == 2

thetaguess = [ 0.2408
               0.8294
               0.2592
               2.4926
               0.8463];
         
elseif answcase == 3
thetaguess = [ 0.1344
               0.6899
               0.6159
               0.5624
               0.9322];
end
% rephrasing the guesses about correlation on the logit scale
% this is done to help the optimizing routine not run into trouble
thetaguess(2) = log(thetaguess(2)/(1-thetaguess(2)));
thetaguess(5) = log(thetaguess(5)/(1-thetaguess(5)));

% computing sum of squares deviations at the initial guess for parameters:
% the function *mom* (method of moments) simulates 500 samples of 250 periods each,
% compoutes relevant moments about these simulated data and compares how far they are from the 
% targeted moments and returns this distance as "sumofsquares"

init = mom(thetaguess,poids,moments);
disp([ 'sumofsquares at initial guess: ' num2str(init)])
   
% Send to optimizing routines using the Simplex algorithm:
% the routine will repeatedly call "mom" with different values for the
% parameter vector and try to minimize "sumofsquares"
options = optimset('display','iter','MaxFunEvals',400,'TolX',1e-2,'TolFun',1e-2,'PlotFcns',@optimplotfval); 
[theta, fval] = fminsearch('mom',thetaguess,options,poids,moments);

disp('print the ending parameter vector in long format:')
format long
disp(theta)
format short

% when finished make sure I got the correct figures by running mom one last
% time
final = mom(theta,poids,moments);


% report the final estimate using the notation from Table 3
theta(2) = 1/(1+exp(-theta(2)));
theta(5) = 1/(1+exp(-theta(5)));

sige_t = 0.01;
sige_p = theta(1)*0.01; % std. deviation of innovations to those shocks
rhoe_t = theta(2); % serial correlation of transitory oil-supply shocks

sigz_t = theta(3)*0.01;
sigz_p = theta(4)*theta(3)*0.01; % std. deviation of innovations to those shocks
rhoz_t = theta(5); 

    disp(' ')
    disp('********************************************************')
    disp('Final Parameter vector, in the notation used in Table 3: ')
    disp([ ' sige_t = ' num2str(sige_t) ])
    disp([ ' sige_p = ' num2str(sige_p) ])
    disp([ ' rhoe_t = ' num2str(rhoe_t) ])
    disp([ ' sigz_t = ' num2str(sigz_t) ])
    disp([ ' sigz_p = ' num2str(sigz_p) ])
    disp([ ' rhoz_t = ' num2str(rhoz_t) ])
  


